SAMSHA NSSATS Data Cleaning

library(dplyr)
library(tidyr)

#could i look at PRIMPAY and HLTHINS columns to investigate moral hazard on the part of treatment centers for patients with certain insurance/payment characteristics (longer stays)?

admission_cols <- c("CASEID", "PREG", "STFIPS",  # state ID 
                    "SERVICES",  # type of treatment / service setting
                    "PSOURCE",  # method of referral to treatment
                    "NOPRIOR", # prior admissions to treatment
                    "SUB1", "ROUTE1", "FREQ1", "FRSTUSE1", # characteristics of their primary substance, method and frequency of use  (could also include second substances)
                    "HLTHINS", "PRIMPAY", # payment details of interest
                    "REASON", # reason for discharge
                    "LOS" #length of stay 
                    )

nssumhs2021 <- read.csv("../Data/NSUMHSS/NSUMHSS_2021_PUF_CSV.csv")
nssumhs2022 <- read.csv("../Data/NSUMHSS/NSUMHSS_2022_PUF_CSV.csv")
nssumhs2023 <- read.csv("../Data/NSUMHSS/NSUMHSS_2023_PUF_CSV.csv")

head(nssumhs2021)

table(nssumhs2021$OWNERSHP) #M and L are missing data 

    1     2     3     4     5     6     L     M 
 7370 11828   735  1109   215   600   479 10035 

table(nssumhs2021$TREATMT_SU) #according to documentation, the TREATMT_SU question asks whether they offer Substance use treatment

    0     1     L     M 
  939 16776  4687  9969 

#we can filter to where the answer is "1" for yes (not 0 for no, L for logical skip, nor M for missing) 

clean_nssumhs <- function(df){
  clean_df <- df %>% 
    #filter to where TREATMT_SU is 1 , meaning they provide substance use treatment
  filter(TREATMT_SU == 1) %>%
    #rename LOCATIONSTATE to STATE to match NSSATS
  rename(STATE = LOCATIONSTATE) %>%
    #replace M and L in ownership column with missing (na)
  #mutate(OWNERSHP = ifelse(is.na(OWNERSHP), NA, OWNERSHP)) %>%
    #alternatively can filter out missing ownership columns
    filter((OWNERSHP != "M") & (OWNERSHP != "L")) %>%
    #cast ownership to integer
    mutate(OWNERSHP = as.integer(OWNERSHP))
  
  return(clean_df)
}
clean_nssumhs2021 <- clean_nssumhs(nssumhs2021)
clean_nssumhs2022 <- clean_nssumhs(nssumhs2022)
clean_nssumhs2023 <- clean_nssumhs(nssumhs2023)

#nssats2009 <- read.csv("N-SSATS-2009-DS0001-bndl-data-tsv/N-SSATS-2009-DS0001-data/N-SSATS-2009-DS0001-data-excel.csv")
nssats2014 <- read.csv("../Data/NSSATS/N-SSATS-2014-DS0001-bndl-data-tsv/N-SSATS-2014-DS0001-data/N-SSATS-2014-DS0001-data-excel.csv")
nssats2015 <- read.csv("../Data/NSSATS/N-SSATS-2015-DS0001-bndl-data-tsv/N-SSATS-2015-DS0001-data/N-SSATS-2015-DS0001-data-excel.csv")
nssats2016 <- read.csv("../Data/NSSATS/NSSATSPUF_2016.csv")
nssats2017 <- read.csv("../Data/NSSATS/NSSATS_PUF_2017_CSV.csv")
nssats2018 <- read.csv("../Data/NSSATS/NSSATS_PUF_2018_CSV.csv")
nssats2019 <- read.csv("../Data/NSSATS/NSSATS_PUF_2019_CSV.csv")
nssats2020 <- read.csv("../Data/NSSATS/NSSATS_PUF_2020_CSV.csv")
head(nssats2015)

Data cleaning on all dfs

profit_types <- c("Forprofit", "Nonprofit", "Government", "Government", "Government", "Government")

# Function to append a number to column names (from ChatGPT)
append_number <- function(name, num) {
  paste(name, num, sep = "_")
}

#year_dfs <- list(nssats2014, nssats2015, nssats2016, nssats2017, nssats2018, nssats2019)#, nssats2020)
year_dfs <- list(nssats2014, nssats2015, nssats2016, nssats2017, nssats2018, nssats2019, nssats2020, 
                 clean_nssumhs2021, clean_nssumhs2022, clean_nssumhs2023)

#year_nums <- seq(2014, 2019)
year_nums <- seq(2014, 2023)

results <- list()
index <- 1
for (year_df in year_dfs)
{
  this_df <- year_df %>% 
    mutate(type = profit_types[OWNERSHP]) %>%
    filter(! STATE %in% c("VI", "GU", "PW", "PR", "FM")) %>% #filter out US territories
    group_by(STATE, type) %>% 
    summarize(count = n()) %>%
    ungroup() %>%
    group_by(STATE) %>%
    mutate(percent = count / sum(count)) %>% 
    ungroup() %>%
    pivot_wider(id_cols = STATE, names_from = type, values_from = c("count", "percent")) %>% #pivot long rows to columns
    rename_with(~ ifelse(. == "STATE", ., append_number(., year_nums[index])), everything())  #add the year onto the columns to track %>%
  

  # add the dataframe to a list 
  results[[index]] <- this_df
  index <- index + 1
  
}

`summarise()` has grouped output by 'STATE'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE'. You can override using the `.groups` argument.

Merge the dataframes together

# profit_states_1419 <- results[[1]]
# 
# for (index in seq(2, length(results)))
# {
#   profit_states_1419 <- left_join(profit_states_1419, results[[index]])
# }
# 
# #add columns that calculate the change from 2015 - 2019 
# profit_states_1419$change_Forprofit_1519 <- profit_states_1419$count_Forprofit_2019 - profit_states_1419$count_Forprofit_2015
# profit_states_1419$change_Nonprofit_1519 <- profit_states_1419$count_Nonprofit_2019 - profit_states_1419$count_Nonprofit_2015
# profit_states_1419$change_Government_1519 <- profit_states_1419$count_Government_2019 - profit_states_1419$count_Government_2015

profit_states_1423 <- results[[1]]

for (index in seq(2, length(results)))
{
  profit_states_1423 <- left_join(profit_states_1423, results[[index]])
}

Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`

#add columns that calculate the change from 2015 - 2019 
profit_states_1423$change_Forprofit_1523 <- profit_states_1423$count_Forprofit_2023 - profit_states_1423$count_Forprofit_2015
profit_states_1423$change_Nonprofit_1523 <- profit_states_1423$count_Nonprofit_2023 - profit_states_1423$count_Nonprofit_2015
profit_states_1423$change_Government_1523 <- profit_states_1423$count_Government_2023 - profit_states_1423$count_Government_2015

#merge state abbreviations to state names
state_names <- tolower(c('ALABAMA','ALASKA','AMERICAN SAMOA','ARIZONA','ARKANSAS','CALIFORNIA','COLORADO','CONNECTICUT','DELAWARE','DISTRICT OF COLUMBIA','FLORIDA','GEORGIA','GUAM','HAWAII','IDAHO','ILLINOIS','INDIANA','IOWA','KANSAS','KENTUCKY','LOUISIANA','MAINE','MARYLAND','MASSACHUSETTS','MICHIGAN','MINNESOTA','MISSISSIPPI','MISSOURI','MONTANA','NEBRASKA','NEVADA','NEW HAMPSHIRE','NEW JERSEY','NEW MEXICO','NEW YORK','NORTH CAROLINA','NORTH DAKOTA','NORTHERN MARIANA IS','OHIO','OKLAHOMA','OREGON','PENNSYLVANIA','PUERTO RICO','RHODE ISLAND','SOUTH CAROLINA','SOUTH DAKOTA','TENNESSEE','TEXAS','UTAH','VERMONT','VIRGINIA','VIRGIN ISLANDS','WASHINGTON','WEST VIRGINIA','WISCONSIN','WYOMING'))
state_abbrevs <- c('AL','AK','AS','AZ','AR','CA','CO','CT','DE','DC','FL','GA','GU','HI','ID','IL','IN','IA','KS','KY','LA','ME','MD','MA','MI','MN','MS','MO','MT','NE','NV','NH','NJ','NM','NY','NC','ND','MP','OH','OK','OR','PA','PR','RI','SC','SD','TN','TX','UT','VT','VA','VI','WA','WV','WI','WY')
state_abbrevs_df <- data.frame(state_names, state_abbrevs)

#we may want to switch to using the Freds population data so it is consistent and we can include from 2013-2024 MERGING DATA TOGETHER TO INCLUDE COVARIATES:

state_colnames <- c("State_Name", "Y2012", "Y2013", "Y2014", "Y2015", "Y2016", "Y2017", "Y2018", "Y2019", "Y2020")

#statepop_long <- read.csv("../Data/state_health_spending/US_POPULATION20.CSV") %>% 
#  dplyr::select(state_colnames) %>%
#  pivot_longer(cols = state_colnames[-1]) %>% 
#  filter(State_Name != "") %>%
#  mutate(year = gsub("Y", "", name)) %>%
#  rename(pop_thousands = value) %>%
#  dplyr::select(State_Name, year, pop_thousands)

#statepop_long

state_healthspending_long <- read.csv("../Data/state_health_spending/US_PER_CAPITA20.CSV") %>%
  dplyr::select(state_colnames, "Item") %>%
  pivot_longer(cols = state_colnames[-1]) %>% 
  #filter to Personal Health Care, which is total spending 
  filter(State_Name != "", grepl("Personal Health Care", Item)) %>%
  mutate(year = gsub("Y", "", name)) %>%
  rename(percapita_health_spending = value, spending_type = Item) %>%
  dplyr::select(State_Name, year, percapita_health_spending, spending_type)

state_healthspending_long 

question - will medicaid enrollment be correlated with population? should we divide by that?

state_medicaidenroll_long <- read.csv("../Data/state_health_spending/MEDICAID_ENROLLMENT20.CSV") %>%
   dplyr::select(state_colnames) %>%
   pivot_longer(cols = state_colnames[-1]) %>% 
   filter(State_Name != "") %>%
   mutate(year = gsub("Y", "", name)) %>%
   rename(medicaid_enrollment = value) %>%
   dplyr::select(State_Name, year, medicaid_enrollment)

state_medicaidenroll_long

# clean_nsduh_files <- function(folder_path, measure_name)
# {
#   #create columns to save the confidence intervals in 
#   measure_upper <- paste(measure_name, "_ciupper", sep = "")
#   measure_lower <- paste(measure_name, "_cilower", sep = "")
#   
#   #loop through all data in the folder, reformat, and save it
#   year_data <- list()
#   nsduh_files <- list.files(path = folder_path)
#   print(nsduh_files)
#   index <- 1
#   
#   #loop through all the files in the folder 
#   for (file in nsduh_files)
#   {
#       this_data <- read.csv(paste(folder_path, file, sep = "")) %>%
#           #round the year range (i.e. 2013-2014) up, so we are saving it as the year_end
#           separate(year_pair, into = c(NA, "year_end")) %>%
#           mutate(year_end = paste("20", year_end, sep = "")) %>%
#           rename(!!measure_name := estimate,
#                  !!measure_lower := ci_lower,
#                  !!measure_upper := ci_upper) %>% 
#           dplyr::select(-c(outcome, age_group))
#       
#       #add it to a list
#       year_data[[index]] <- this_data
#       index <- index + 1
#   }
# 
#   #append the data from each file into one dataframe and return
#   return(bind_rows(year_data))
# }
# 
# 
# sud_pastyear <- clean_nsduh_files("Data/NSDUH/NSDUH_SUD/", "sud_pastyear")
# needing_notreceiving <- clean_nsduh_files("Data/NSDUH/NSDUH_NeedingNotReceiving/", "needing_tx")
# bingealcohol_pastmonth <- clean_nsduh_files("Data/NSDUH/NSDUH_BingeAlcohol/", "bingealcohol_pastmonth")
# drug_pastmonth <- clean_nsduh_files("Data/NSDUH/NSDUH_IllicitDrugUse/", "drug_pastmonth")
# 
# sud_demand_df <- sud_pastyear %>% 
#   left_join(needing_notreceiving, by= c("state", "year_end")) %>%
#   left_join(bingealcohol_pastmonth, by= c("state", "year_end")) %>%
#   left_join(drug_pastmonth, by = c("state", "year_end"))
# 
# sud_demand_df

Read all NSDUH data from here: https://www.samhsa.gov/data/nsduh/state-reports

library(haven)
nsduh_99_19 <- haven::read_sas("../Data/NSDUH/NSDUH_99_19_state_saes_final.sas7bdat")


nsduh_21_23 <- haven::read_sas("../Data/NSDUH/state_saes_final_2122_22_2223.sas7bdat")

#hitting a hiccup with 2021-2023 NSDUH data - documentation says that the previous years cannot be compared to 2021-2023
# https://www.samhsa.gov/data/sites/default/files/reports/rpt56198/2023-nsduh-puf-data-users-guide.pdf
#nsduh_2020 <- read.table("../Data/NSDUH/NSDUH_2020_Tab.txt", header = TRUE, sep = '\t')
#nsduh_2021 <- read.table("../Data/NSDUH/NSDUH_2021_Tab.txt", header = TRUE, sep = '\t')
#nsduh_2022 <- read.table("../Data/NSDUH/NSDUH_2022_Tab.txt", header = TRUE, sep = '\t')
#nsduh_2023 <- read.table("../Data/NSDUH/NSDUH_2023_Tab.txt", header = TRUE, sep = '\t')

looking for variables which can measure demand but are consistently available for a longer period of time - alcohol use disorder in the past year, (ABODALC) cocaine use in the past year, (COCYR) heroin use in the past year, (HERYR) marijuana use in the past month (MRJMON) - month chosen over year because it is more frequent and more likely represents high use / use disorder

(probably won’t include tobacco use because it isnt as often treated in a specialty facility)

clean_nsduh_data <- function(nsduh_data){
 #filter the data to only state reports, and years post-2015 
nsduh_data_long <- nsduh_data %>% 
  separate(pyearnm, into = c("start_year", "end_year"), "-") %>%
  filter(area == 2) %>% 
  filter(outcome %in% nsduh_outcomes) %>%
  filter(agegrp  == 4) %>% #filter to people above 18 (adults)
  mutate(stname_lower = tolower(stname)) %>%
  dplyr::select(outname, outcome, start_year, end_year, agegrp, stname,stname_lower, BSAE, low_sae, up_sae)
                #pop, est_total, low_total, up_total) #i was previously using est_total but this is the average number of people per state which will be dependent on population

nsduh_data_wide <- nsduh_data_long %>%
  select(stname, start_year, end_year, outcome, BSAE) %>%
  pivot_wider(id_cols = c(stname, end_year), names_from = outcome, values_from = BSAE, unused_fn = first)

  return(nsduh_data_wide) 
}
#interested in 4 substance use outcomes
nsduh_outcomes <- c("BNGDRK", "ILLEMMON", "TXNPILAL", "UDPYILAL", "HERYR", "ABODALC", "UDYR5ALC", "COCYR", "HERYR", "MRJMON")

nsduh_data_wide <- clean_nsduh_data(nsduh_99_19)
nsduh_data_wide_21 <- clean_nsduh_data(nsduh_21_23)

Warning: Expected 2 pieces. Missing pieces filled with `NA` in 1064 rows [8345, 8346, 8347, 8348, 8349, 8350, 8351, 8352, 8353, 8354, 8355, 8356, 8357, 8358, 8359, 8360, 8361, 8362, 8363, 8364, ...].

nsduh_data_wide

nsduh_data_wide_21

NA

#this function assumes that statepop_long, state_healthspending_long, and state_medicaidenroll_long are read into memory already
merge_control_vars <- function(nsduh_data, alc_colname)
{
  
  nsduh_data$stname <- trimws(nsduh_data$stname)
  
  state_control_long <- state_healthspending_long %>%
    #statepop_long %>% 
  #inner_join(state_healthspending_long, by = c("State_Name"= "State_Name", "year" = "year")) %>%
  inner_join(state_medicaidenroll_long, by = c("State_Name"= "State_Name", "year" = "year")) %>%
  #mutate(medicaid_enroll_percapita = medicaid_enrollment / pop_thousands) %>%
  
    
    #join "demand" measures from NSDUH data
  #left_join(nsduh_data, by = c("State_Name" = "stname", "year" = "end_year")) %>% 
  #changing to right join on 6/12 because we don't have control data for 2021-2023 yet
  right_join(nsduh_data, by = c("State_Name" = "stname", "year" = "end_year")) %>% 
  #join state abbreviations onto dataframe
  mutate(state_name_lower = tolower(State_Name)) %>%
  left_join(state_abbrevs_df, by = c("state_name_lower" = "state_names")) %>%
  dplyr::select(-c(spending_type, state_name_lower)) %>%
  rename("binge_drinking" = BNGDRK, "illicit_druguse" = ILLEMMON, 
         #"needing_tx"= TXNPILAL, 
         #"sud" = UDPYILAL,
         "aud_yr" = {{alc_colname}}, 
         "heroin_yr"=HERYR, "cocaine_yr"=COCYR, "marijuana_month" = MRJMON)

  return(state_control_long)
}

state_control_long14 <- merge_control_vars(nsduh_data_wide, "ABODALC")
state_control_long21 <- merge_control_vars(nsduh_data_wide_21, "UDYR5ALC")

state_control_long <- bind_rows(state_control_long14, state_control_long21)
state_control_long <- state_control_long %>% filter(start_year > 2012)
state_control_long

Scale the variables and lag the need variables

#scale to health spending in thousands of dollars
state_control_long$percapita_health_spending_thousands <- state_control_long$percapita_health_spending / 1000
#scale to population in millions 
#state_control_long$pop_millions <- state_control_long$pop_thousands / 1000

#change the demand variables to percents
state_control_long$binge_drinking <- state_control_long$binge_drinking * 100
state_control_long$illicit_druguse <- state_control_long$illicit_druguse * 100
#state_control_long$sud <- state_control_long$sud * 100 
#state_control_long$needing_tx <- state_control_long$needing_tx * 100
state_control_long$aud_yr <- state_control_long$aud_yr * 100
state_control_long$heroin_yr <- state_control_long$heroin_yr * 100
state_control_long$cocaine_yr <- state_control_long$cocaine_yr * 100
state_control_long$marijuana_month <- state_control_long$marijuana_month * 100

#lag the 4 demand variables , and change them to percentages 
state_control_long <- state_control_long %>%
                    arrange(State_Name, year) %>%  # Ensure data is ordered by State_Abbrev and year
                    group_by(State_Name) %>%
                    mutate(prev2years_bingedrinking = dplyr::lag(binge_drinking),
                           prev2years_druguse = dplyr::lag(illicit_druguse),
                           #prev2years_sud = dplyr::lag(sud),
                           #prev2years_needingtx = dplyr::lag(needing_tx),
                           prev2years_aud = dplyr::lag(aud_yr),
                           prev2years_marijuana = dplyr::lag(marijuana_month),
                           druguse_calc = (cocaine_yr + heroin_yr)/2, #combine the 2 drug use metrics - average 
                           #BWO 5/22- changing the above metric to divide by 2 for the prevalence of drug use per 100
                           prev2years_drugusecalc = dplyr::lag(druguse_calc),
                           prev2years_cocaine = dplyr::lag(cocaine_yr),
                           prev2years_heroin = dplyr::lag(heroin_yr))

#pivot data to long format
#state_final_df <- profit_states_1419 %>%
state_final_df <- profit_states_1423 %>%
  dplyr::select(STATE, contains("count")) %>%
  pivot_longer(cols = contains("count")) %>%
  separate(name, into = c(NA, "profit_type", "year"), "_") %>%
  rename(State_Abbrev = STATE, count_tx = value) %>%
  mutate(State_Abbrev = as.character(State_Abbrev), year = as.character(year)) %>%
#join to control variables
  left_join(state_control_long, by = c("State_Abbrev" = "state_abbrevs", "year" = "year")) #%>%
  #filter(year >= 2016 & year <= 2020) %>%
  #mutate(tx_percapita = count_tx / (pop_thousands*1000))

head(state_final_df)

Merging following variables: * race * gender * ethnicity

#read demographic information from 1 year from the WONDER CDC dataset
read_clean_demographic <- function(filepath, year_param)
{
  demographic_states <- read.table(filepath, sep = "\t", 
           fill = TRUE,
           header = TRUE) %>%
    #after 2020 it is a "Sex" variable instead of "Gender"
    rename(any_of(c(Gender = "Sex"))) %>%
  filter(!is.na(Population),#filter out unnecessary text at the bottom 
         Gender != "", #filter to the most granular level of each combination of gender/Race/Ethnicity and then we can aggregate up
         Race != "", 
         Ethnicity != "") 

gender_df <- demographic_states %>%  
  arrange(States) %>% 
  group_by(States, Gender) %>% 
  summarize(count = sum(Population)) %>%
  #change the count per gender into a percent
  mutate(Total = sum(across(where(is.numeric))),
         Percent = (count / Total) * 100)  %>% 
   ungroup() %>%
  pivot_wider(id_cols = States, names_from = Gender, values_from = Percent, names_glue = "Gender_{Gender}") 

race_df <- demographic_states %>%  
  arrange(States) %>% 
  group_by(States, Race) %>% 
  summarize(count = sum(Population)) %>%
  #change the count per Race into a percent
  mutate(Total = sum(across(where(is.numeric))),
        #change the percent from a decimal (i.e. 10% = 0.1) to a percent (i.e 10% = 10) to correspond with other variables' format in df
         Percent = (count / Total) * 100) %>% 
   ungroup() %>%
  #replace all spaces in the Race column
  mutate(Race = gsub(" ", "", Race)) %>% 
  #Pivot the race into columns
  pivot_wider(id_cols = States, names_from = Race, values_from = Percent, names_glue = "Race_{Race}") 

ethnicity_df <- demographic_states %>%  
  arrange(States) %>% 
  group_by(States, Ethnicity) %>% 
  summarize(count = sum(Population)) %>%
  #change the count per Ethnicity into a percent
  mutate(Total = sum(across(where(is.numeric))),
         Percent = (count / Total) * 100)  %>% 
   ungroup() %>%
  #replace all spaces in the Ethnicity column
  mutate(Ethnicity = gsub(" ", "", Ethnicity)) %>%
  #pivot the ethnicity into columns
  pivot_wider(id_cols = States, names_from = Ethnicity, values_from = Percent, names_glue = "Ethnicity_{Ethnicity}") 

demographics_df <- gender_df %>% 
  merge(ethnicity_df, by = "States") %>%
  merge(race_df, by = "States") %>%
  mutate(year = year_param)

return(demographics_df)
}


demographics_2013 <- read_clean_demographic("../Data/StateYear_DemographicInfo/Single-Race Population Estimates 2010-2020 by State and Single-Year Age_2013.txt", 2013)

`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.

demographics_2014 <- read_clean_demographic("../Data/StateYear_DemographicInfo/Single-Race Population Estimates 2010-2020 by State and Single-Year Age_2013.txt", 2014)

`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.

demographics_2015 <- read_clean_demographic("../Data/StateYear_DemographicInfo/Single-Race Population Estimates 2010-2020 by State and Single-Year Age_2013.txt", 2015)

`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.

demographics_2016 <- read_clean_demographic("../Data/StateYear_DemographicInfo/Single-Race Population Estimates 2010-2020 by State and Single-Year Age_2013.txt", 2016)

`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.

demographics_2017 <- read_clean_demographic("../Data/StateYear_DemographicInfo/Single-Race Population Estimates 2010-2020 by State and Single-Year Age_2013.txt", 2017)

`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.

demographics_2018 <- read_clean_demographic("../Data/StateYear_DemographicInfo/Single-Race Population Estimates 2010-2020 by State and Single-Year Age_2013.txt", 2018)

`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.

demographics_2019 <- read_clean_demographic("../Data/StateYear_DemographicInfo/Single-Race Population Estimates 2010-2020 by State and Single-Year Age_2019.txt", 2019)

`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.

demographics_2020 <- read_clean_demographic("../Data/StateYear_DemographicInfo/Single-Race Population Estimates 2020-2022 by State and Single-Year Age_2020.xls", 2020)

`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.

demographics_2021 <- read_clean_demographic("../Data/StateYear_DemographicInfo/Single-Race Population Estimates 2020-2022 by State and Single-Year Age_2021.xls", 2021)

`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.

demographics_2022 <- read_clean_demographic("../Data/StateYear_DemographicInfo/Single-Race Population Estimates 2020-2022 by State and Single-Year Age_2022.xls", 2022)

`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.`summarise()` has grouped output by 'States'. You can override using the `.groups` argument.

state_demographics <- bind_rows(demographics_2013,
          demographics_2014,
          demographics_2015,
          demographics_2016,
          demographics_2017,
          demographics_2018, 
          demographics_2019,
          demographics_2020,
          demographics_2021,
          demographics_2022) %>%
  mutate(year = as.character(year)) #cast the year column to character

#delete unnecessary dfs 
rm(demographics_2013, demographics_2014, demographics_2015, demographics_2016, demographics_2017, demographics_2018, demographics_2019, demographics_2020,demographics_2021, demographics_2022)

Merging insurance variable * insurance

read_clean_insurance <- function(path, year_param)
  {
  insurance_states <- read.csv(path,
                               skip=2,
                               header = TRUE) %>%
  filter(!is.na(Total)) %>%
  select(-Footnotes) %>%
  mutate(year = year_param) %>%
  rename_with(~paste0("Insurance_", .x), everything()) %>%
    #remove the "Insurance" from year and Location
  rename(year = Insurance_year, Location = Insurance_Location)
  
  #some of the columns show "<0.01" - calculate them by subtracting from the other values 
  insurance_states <- insurance_states %>%
    mutate(Insurance_Military = na_if(as.character(Insurance_Military), "<.01")) %>% #change any "<.01" text to na
    mutate(Insurance_Military = coalesce(as.numeric(Insurance_Military), 1 - (Insurance_Employer + Insurance_Non.Group + Insurance_Medicaid + Insurance_Medicare + Insurance_Uninsured))) #calculate the remaining percentage
  
  #change the decimals to percents 
  insurance_states <- insurance_states %>% 
                      mutate(year = as.character(year)) %>% # ensure the year is a character column
                      mutate_if(is.numeric, ~. * 100)
  
  return(insurance_states)
}

insurance_2013 <- read_clean_insurance("../Data/StateYear_Insurance/raw_data2013.csv", 2013)
insurance_2014 <- read_clean_insurance("../Data/StateYear_Insurance/raw_data2014.csv", 2014)
insurance_2015 <- read_clean_insurance("../Data/StateYear_Insurance/raw_data2015.csv", 2015)
insurance_2016 <- read_clean_insurance("../Data/StateYear_Insurance/raw_data2016.csv", 2016)
insurance_2017 <- read_clean_insurance("../Data/StateYear_Insurance/raw_data2017.csv", 2017)
insurance_2018 <- read_clean_insurance("../Data/StateYear_Insurance/raw_data2018.csv", 2018)
insurance_2019 <- read_clean_insurance("../Data/StateYear_Insurance/raw_data2019.csv", 2019)
#insurance_2020 <- read_clean_insurance("../Data/StateYear_Insurance/raw_data2020.csv", 2020)
insurance_2021 <- read_clean_insurance("../Data/StateYear_Insurance/raw_data2021.csv", 2021)

Warning: There was 1 warning in `mutate()`.
ℹ In argument: `Insurance_Military = coalesce(...)`.
Caused by warning in `list2()`:
! NAs introduced by coercion

insurance_2022 <- read_clean_insurance("../Data/StateYear_Insurance/raw_data2022.csv", 2022)

state_insurance <- bind_rows(insurance_2013, insurance_2014, insurance_2015, insurance_2016, insurance_2017, insurance_2018, insurance_2019, 
                             #insurance_2020, 
                             insurance_2021, insurance_2022) %>%
  mutate(year = as.character(year))

#delete unnecessary dfs
rm(insurance_2013, insurance_2014, insurance_2015, insurance_2016, insurance_2017, insurance_2018, insurance_2019, #insurance_2020, 
   insurance_2021, insurance_2022)

Merging poverty and unemployment variables

library(janitor) #use to overwrite row names 
library(tibble)
library(readxl)

freds_name_mapping <- read_excel("../Data/StateYear_Income_Unemployed/State_Population_Unemployment_Income.xlsx", sheet=2) 
#create a named vector of the values from the name mapping
state_names <- freds_name_mapping$State_Parsed
names(state_names) <- freds_name_mapping$`Series ID:`

value_names <- freds_name_mapping$`Title:`
names(value_names) <- freds_name_mapping$`Series ID:`

freds_df_unformatted <- read_excel("../Data/StateYear_Income_Unemployed/State_Population_Unemployment_Income.xlsx", sheet=3) 

freds_df <- freds_df_unformatted %>%
  #switch the rows and the columns so the values are in the rows and columns are years
  t() %>%  
  #cast to dataframe
  as.data.frame() %>% 
  #move the first row to be the header
  janitor::row_to_names(row_number = 1) %>% 
  #rename the row names to "freds_colname"
  tibble::rownames_to_column("freds_colname") %>% 
  mutate(val = value_names[freds_colname]) %>%
  #mark which variables we are looking at 
  separate(val, into = c("variable_measure", "state_name"), sep = " in ") %>% 
  #remove "the" from DC name
  mutate(state_name = ifelse(state_name == "the District of Columbia", "District of Columbia", state_name)) %>% 
  #shorten variable measure names 
  mutate(variable_measure = ifelse(variable_measure == "Resident Population", "Population", 
                                   ifelse(variable_measure == "Unemployment Rate", "Unemployment", 
                                          ifelse(variable_measure == "Real Median Household Income", "Income", "MISSING")))) %>%
  #remove month/date from year names - unnecessary
  rename("2013" = "2013-01-01",
         "2014" = "2014-01-01",
         "2015" = "2015-01-01",
         "2016" = "2016-01-01",
         "2017" = "2017-01-01",
         "2018" = "2018-01-01",
         "2019" = "2019-01-01",
         "2020" = "2020-01-01",
         "2021" = "2021-01-01",
         "2022" = "2022-01-01",
         "2023" = "2023-01-01") %>% 
  #change the years from the columns to a row indicator - data is now in long format with one row per state-year-measurement
  pivot_longer(cols = c("2013", "2014", "2015", "2016", "2017", "2018", "2020", "2019", "2022", "2023")) %>%
  rename("year" = name) %>%
  #pivot the individual measures into columns so we have only one row per state-year
  pivot_wider(id_cols = c(state_name, year), values_from = value, names_from = variable_measure) %>% 
  #cast numeric variables to numeric
  mutate(year = as.character(year), Population = as.numeric(Population), Unemployment = as.numeric(Unemployment), Income = as.numeric(Income))

state_final_df <- state_final_df %>%
  #join demographic information
  left_join(state_demographics, by = c("State_Name"="States", "year"="year")) %>%
  #join insurance information
  left_join(state_insurance, by = c("State_Name"="Location", "year"="year")) %>%
  #join the population, unemployment, and income information
  left_join(freds_df, by=c("State_Name"="state_name", "year"="year"))

Error: object 'state_final_df' not found

These are available starting at 2013-2014: alcohol use disorder in the past year, cocaine use in the past year, heroin use in the past year, marijuana use in the past month (probably won’t include tobacco use because it isnt as often treated in a specialty facility)

# state_forprofit <- state_final_df %>% 
#                     arrange(State_Abbrev, year) %>%  # Ensure data is ordered by State_Abbrev and year
#                     group_by(State_Abbrev) %>%
#                     filter(profit_type == "Forprofit", !is.na(prev2years_drugusecalc)) #2014 does not have this data
# 
# state_nonprofit <- state_final_df %>% 
#                     arrange(State_Abbrev, year) %>%  # Ensure data is ordered by State_Abbrev and year
#                     group_by(State_Abbrev) %>%
#                     filter(profit_type == "Nonprofit", !is.na(prev2years_drugusecalc))
# 
# state_govt <- state_final_df %>% 
#                     arrange(State_Abbrev, year) %>%  # Ensure data is ordered by State_Abbrev and year
#                     group_by(State_Abbrev) %>%
#                     filter(profit_type == "Government", !is.na(prev2years_drugusecalc))
# 
# 
# state_totals <- state_final_df %>%
#                   filter(!is.na(prev2years_drugusecalc)) %>%  #filter out rows we don't want first 
#                   group_by(State_Abbrev, year) %>% 
#                   mutate(count_tx_total = sum(count_tx, na.rm = TRUE)) %>% #calculate the total number of treatment centers per year 
#                   select(-c(count_tx, profit_type)) %>%  #remove the columns that only applied to specific treatment center counts 
#                   ungroup() %>%
#                   rename(count_tx = count_tx_total) %>%
#                   distinct(State_Abbrev, year, count_tx, .keep_all = TRUE)# get rid of the duplicate rows of control variables for each profit status
# 
# write.csv(state_forprofit, "Data/CleanData/Forprofit_finaldata.csv", row.names = FALSE)
# write.csv(state_nonprofit, "Data/CleanData/Nonprofit_finaldata.csv", row.names = FALSE)
# write.csv(state_govt, "Data/CleanData/Government_finaldata.csv", row.names = FALSE)
# write.csv(state_totals, "Data/CleanData/Totals_finaldata.csv", row.names = FALSE)

state_final_df <- state_final_df %>%
  pivot_wider(id_cols = c(State_Abbrev, year, "State_Name", "pop_millions", "percapita_health_spending_thousands", "medicaid_enroll_percapita", "prev2years_drugusecalc", "prev2years_aud", "prev2years_marijuana", "Gender_Male", "Ethnicity_HispanicorLatino", "Race_White", "Insurance_Uninsured", "Unemployment", "Income"), #keep all variables we will need for the model
              names_from = profit_type, values_from = count_tx) %>% 
  #calculate the sum of all profit types
  rowwise() %>%
          mutate(count_tx_total = sum(Forprofit, Nonprofit, Government, na.rm = TRUE)) %>%
  #filter only to data that will allow the lagging 
  filter(year >= 2015) %>%
  #replace the nan in Government with 0
  mutate(Government = ifelse(is.na(Government), 0, Government))
  

drugoverdose_df <- read.csv("../Data/CDC_DrugOverdoseDeaths.csv")

#format the data 
drugoverdose_df <- drugoverdose_df %>%
  arrange(STATE, YEAR) %>%
  mutate(YEAR = as.character(YEAR),
         overdose_rate_lag = lag(RATE)) %>%
  rename(overdose_rate = RATE) %>%
  select(STATE, YEAR, overdose_rate, overdose_rate_lag)


#merge the overdose data onto the final df
state_final_df <- state_final_df %>% 
  left_join(drugoverdose_df, by = c('year'='YEAR', 'State_Abbrev'='STATE'))

write.csv(state_final_df, "../Data/CleanData/Allprofit_finaldata_6.12.25.csv", row.names = FALSE)

As of 6/12/25 - creating a new dataset that separates out the OTPs:

nssats2014 %>% filter(OTP==1)

OTP_results <- list()
OTP_index <- 1
for (year_df in year_dfs)
{
  this_df <- year_df %>% 
    mutate(type = profit_types[OWNERSHP],
           OTP_binary = ifelse(OTP == 1, "OTP", "Not_OTP")) %>%
    filter(! STATE %in% c("VI", "GU", "PW", "PR", "FM")) %>% #filter out US territories
    group_by(STATE, OTP_binary, type) %>% 
    summarize(count = n()) %>%
    mutate(OTP_profit_group = paste(type,OTP_binary, sep = "_")) %>%
    ungroup() %>%
    group_by(STATE) %>%
    mutate(percent = count / sum(count)) %>% 
    ungroup() %>%
    pivot_wider(id_cols = STATE, names_from = OTP_profit_group, values_from = c("count", "percent")) %>% #pivot long rows to columns
    rename_with(~ ifelse(. == "STATE", ., append_number(., year_nums[OTP_index])), everything())  #add the year onto the columns to track %>%
  

  # add the dataframe to a list 
  OTP_results[[OTP_index]] <- this_df
  OTP_index <- OTP_index + 1
  
}

`summarise()` has grouped output by 'STATE', 'OTP_binary'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE', 'OTP_binary'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE', 'OTP_binary'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE', 'OTP_binary'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE', 'OTP_binary'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE', 'OTP_binary'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE', 'OTP_binary'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE', 'OTP_binary'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE', 'OTP_binary'. You can override using the `.groups` argument.`summarise()` has grouped output by 'STATE', 'OTP_binary'. You can override using the `.groups` argument.

profit_states_OTP_1423 <- OTP_results[[1]]

for (OTP_index in seq(2, length(OTP_results)))
{
  profit_states_OTP_1423 <- left_join(profit_states_OTP_1423, OTP_results[[OTP_index]])
}

Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`Joining with `by = join_by(STATE)`

profit_states_OTP_1423